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Abstract. We study localised activity patterns in neural field equations posed on the Euclidean 

plane; such models are commonly used to describe the coarse-grained activity of large ensembles of 

cortical neurons in a spatially continuous way. We employ matrix-free Newton-Krylov solvers and 

perform numerical continuation of localised patterns directly on the integral form of the equation. 

This opens up the possibility to study systems whose synaptic kernel does not lead to an equivalent 

PDE formulation. We present a numerical bifurcation study of localised states and show that the 

proposed models support patterns of activity with varying spatial extent through the mechanism of 

^^^ homoclinic snaking. The regular organisation of these patterns is due to spatial interactions at a 

^^ specific scale associated with the separation of excitation peaks in the chosen connectivity function. 

^^ The results presented form a basis for the general study of localised cortical activity with inputs and, 

^O more specifically, for investigating the localised spread of orientation selective activity that has been 

observed in the primary visual cortex with local visual input. 

^^ Key words, numerical continuation, neural fields, bifurcation, snaking, pattern formation 
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1. Introduction. Oneof the most challenging research questions in neuroscience 

ry^ is understanding the relationship between spatially-structured cortical states and the 

underlying neural circuitry that supports them. A popular approach for analysing 
coarse-grained activity of large ensembles of neurons in the cortex is to model corti- 

F^ cal space as a continuum. Since the pioneering work of Wilson and Cowan [58, 59] 

C^ and Amari [1,2], continuous neural field models have become a popular and effective 

g tool in neuroscience. In such models, the large-scale activity of spatially-extended 

I— I networks of neurons is described in terms of nonlinear integro-differential equations, 

whose associated integral kernels represent the spatial distribution of neuronal synap- 
tic connections. The canonical Wilson- Cowan- Amari neural field equation [2,59] 
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^ ^u(y, t) = -u{r, t)^ j^ w(v, v')S{u(v', t))dm(v') (1.1) 

,/ describes the evolution of the average membrane voltage potential of a neuronal pop- 

^^ ulation ix(r,t), at position r G ^^ on the cortex and at time t. The nonlinear function 

(T^ S represents the neural firing rate, whereas the connectivity function w{y^v') models 

^"H how a population of neurons at position r on the cortex interacts with a population 

^ at position r^ Frequently- used firing rate functions S include the Heaviside step 

•^ function, piecewise-linear functions or smooth sigmoidal functions. Various choices 

rN are also possible for the connectivity function, which is often assumed to be transla- 

^ tion invariant (that is, dependent on the Euclidean distance ||r — r^||) and localised 

in space. The cortical domain Vt is usually a subset of W^^ with (i = 1 or, in more 

realistic models, d = 2. For a recent review on neural fields modeling, we refer to 

Bressloff [8]. 

Unlike spiking neural network models, continuous field models have the advan- 
tage that analytic techniques for partial differential equations (PDEs) can be adapted 
to study the formation of patterns and their dependence upon control parameters. 
Various types of coherent structures have been observed in neural field models, rang- 
ing from spatially and temporally periodic patterns to travelling waves and spiral 
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waves [19,27,42]. Neural field equations have also successfully been used to model a 
wide range of neurobiological phenomena such as visual hallucinations [9,28], mech- 
anisms for short term memory [44] and feature selectivity in the visual cortex [6,35]. 

A common strategy to derive analytical and numerical results for the nonlo- 
cal Eq. (1.1) is to assume translation invariance and exploit the freedom in the choice 
of the connectivity function w: if the Fourier Transform of w is sl rational function, it 
is possible to derive a PDE formulation that is equivalent to the integral model [42]. 
Coherent structures supported by the original model can then be conveniently con- 
structed and analyzed in the PDE framework. Indeed, previous studies have been 
carried out in cases where the synaptic kernel led to an equivalent PDE formula- 
tion [43,44]. To the best of the author's knowledge, there has been no attempt to 
propose efficient path- following methods for general connectivity functions w. This 
paper is motivated by the desire to develop numerical algorithms for Eq. (1.1) without 
relying on an equivalent PDE formulation. More precisely, we discuss how to solve 
Eq. (1.1) when ^^ = R^, 5* is a smooth function of sigmoidal type and w has a generic 
Fourier Transform. 

The main tools for our investigation are time simulation and numerical continu- 
ation. When a PDE model is available, we use standard techniques for both tasks, in 
line with what is typically done in several other works in this field [41-43]. However 
we show that, when the integral of Eq. (1.1) can be written as a convolution, it is 
convenient to employ a fast Fourier transform (FFT) for both time stepping and nu- 
merical continuation. Direct numerical simulations of Eq. (1.1) using FFTs have been 
performed before on full integral models [23,32]. In the present paper, we combine 
FFTs with Newton-Krylov solvers [38,39], thus opening up the possibility to perform 
numerical continuation directly on the integral model. 

We concentrate our effort on the emergence and bifurcation structure of stationary 
localised patterns in planar neural field models of the form Eq. (1.1). Indeed, this type 
of solution is of great interest in models of prefrontal cortex, where localised states 
are believed to characterise working memory [18,33]. Recently, localised regions of 
activity have been observed in the cat primary visual cortex [17] when the animal 
is presented with localised-oriented input. Moreover, some reported drug-induced 
visual hallucinations have also been found to be spatially localised [54] indicating the 
existence of spatially localised regions of activity in the human primary visual cortex. 

Localised states have been observed in a wide variety of nonlinear media [40] . The 
bifurcation structure of localised solutions has been studied extensively in the Swift- 
Hohenberg equation posed on the real line [11,12,16,25,60] and on the plane [3,5, 
45-47,51,52]. In this context, a well-known mechanism for the formation of localised 
states is homodinic snaking: solutions with one or more bumps at the core emerge 
from the trivial homogeneous state and undergo a series of fold bifurcations, giving 
rise to a hierarchy of states with an increasing number of bumps. This scenario seems 
to be a common footprint of localised patterns, extending far beyond the prototypical 
Swift-Hohenberg equation [36,53] and to have also been found in nonlocal equations 
such as neural fields [21,29,43,44]. 

An example of homoclinic snaking is given in Fig. 1.1, where we show a bifurcation 
diagram for the integral model (1.1) posed on the real line with 

w{x, x') := w{\x -x\) = e-^l^-^'l(6sin \x - x'\ + cos(x - x')), (1.2) 

^H = -. ^—^ - T^^ (1-3) 

where b and /i control the decay of the synaptic kernel and the slope of the sigmoidal 
firing rate respectively, while ^ is a threshold value. As fi varies, the trivial steady 
state u = bifurcates at a subcritical Reversible Hopf bifurcation, from which a 
branch of periodic solutions and a branch of localised states originate. Beyond the 
fold Lp on the periodic branch, a stable periodic state, shown in the inset (a), coexists 
with the trivial state. The branch of localised states features solutions with an odd 
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Fig. 1.1. Bifurcation diagram showing snaking behavior for the ID neural field equation us- 
ing Eqs. (1.2)-(1.3) with ^ = 3.5 and h = 0.4^ where stable (unstable) branch segments are rep- 
resented with solid (dashed) lines. Two branches of solutions bifurcate from the trivial states at 
a Reversible Hopf bifurcation EH , namely a branch of periodic solutions (grey) and a branch of 
localised solutions (black). The periodic branch is stable above the fold Lp. The localised branch 
undergoes a series of fold bifurcations giving rise to stable branch segments with increasing numbers 
of bumps. Several examples of solutions are shown in the insets (a)-(d). Stable localised states exist 
for iJ. e [mi,M2]. 



number of bumps and snakes for ja G [/ii,/i2]; in this interval, localised solutions 
with different spatial extent coexist and are stable (see insets (b)-(d)). We note the 
existence of a counterpart even-numbered-bump solution branch along with so-called 
ladder branches connecting the odd and even branches (not shown). For the same 
connectivity function used here, snaking has been shown to occur in terms of the 
parameter b [44]. Elsewhere, Elvin et al. [26] used the Hamiltonian structure of the 
steady states of the model (1.1) with (1.2)-(1.3) and developed numerical techniques 
to find homoclinic orbits of the system. Snaking has also been shown to occur with a 
Mexican-hat connectivity function [21] and for a wizard-hat connectivity function [29]. 
In the latter study, normal form theory for a Reversible Hopf bifurcation was applied 
to prove the existence of localised solutions and a comprehensive parameter study was 
carried out with numerical continuation in terms of two parameters controlling the 
nonlinearity and a third controlling the shape of the connectivity function. 

For the neural field equation posed on the Euclidean plane, various types of 
spatially localised two-dimensional states have been found including radially sym- 
metric solutions [10,30,32,43,44,55,57], rings [23,50], hexagonal patches [30,43] 
along with more complex breathing and travelling states [22,31,50]. When working 
in the Euclidean plane, it is still possible to derive an equivalent PDF for suitable 
choices of the connectivity function w or of its Fourier transform w [43]. Following 
this approach, normal form theory has recently been applied to prove the existence 
of localised solution branches in both the Euclidean plane and on the Hyperbolic 
disk with a wizard-hat connectivity function [30]. When these solutions were path- 
followed using numerical continuation, it was found that the branches do not undergo 
snaking- type behaviour. However, for the connectivity (1.2), snaking was shown to 
occur for branches of radially-symmetric solutions [43]. Furthermore, the existence of 
D6-symmetric and D3-symmetric localised states were found at isolated parameter 
values. 



In Fig. 1.2 we show time simulations of the nonlocal Eq. (1.1) posed on the plane. 
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Fig. 1.2. Time simulations of the integral model (1.1) posed on the plane, with connectivity 
function given by Eq. (I.4) and firing rate function given by Eq. (1.2); in all simulations = 5.6^ 
b = 0.4. (a)-(c): Convergence of a small bump of activity, given by Eq. (A.l), to a spot solution 
with fi = 3.4; time as indicated in panels, (d)-(f): Convergence of a hexagonal pattern, given by 
Eq. (A. 2), to a stable 'D6-symmetric localised state with jjl = 3.2. (g)-(i) With the same initial 
condition as (d), divergence away from a localised state to a periodic state with 11 = 3.4. 



with radially-symmetric connectivity function 



w{v.,v') := w(\\y — y'W) = e 



-6||r 



l(6sin||r-r'||+cos||r-r'||), (1.4) 



and sigmoidal firing rate function given by Eq. (1.3); various combination of initial 
conditions and control parameters lead to three different steady states. We note 
that this models supports localised states, such as the radially-symmetric spot of 
panel (c) or the hexagonal pattern of panel (f). Furthermore, changes in the slope 
of the sigmoidal firing rate affect the stability properties of the solutions, leading to 
other localised states or to domain-covering patterns such as the one in panel (i). In 
the present paper we will focus on localised planar patterns (with various symmetry 
properties) that coexist with the trivial state li = and with fully periodic states, 
similar to the one shown in Fig. 1.1 for the ID case. 

A key result of the present paper is that homoclinic snaking occurs in planar 
neural field models for non-radial patterns and that the choice of the connectivity 
function has a considerable impact on the snaking structure, as well as on the stability 
properties and the selection of the localised states. In two spatial dimensions, periodic 
and localised solution branches bifurcate from the trivial state at a Turing instability 
(as opposed to a Reversible Hopf in ID) and snake irregularly, in a similar fashion to 
what is found for the planar Swift-Hohenberg equation [46]. 

The outline of the paper is as follows. In Sec. 2, we present the different models 
that we study, and then, in Sec. 3, we describe the numerical methods that are used to 
analyze each model. Our homoclinic snaking results of localised states are presented 
in Sec. 4. 
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Fig. 2.1. Firing rate function and connectivity function for the integral model (IM). (a): 
Sigmoidal nonlinearity given by Eq. (2.2) plotted here with = 5.6 and fi = 3; in this formulation 
the effective threshold is 0/fj, (dashed line), (h): Radially- symmetric connectivity kernel given by 
Eq. (2.3) plotted on the Euclidean plane with b = 0.4. (c): Its Fourier transform plotted on the 
{kx iky) -plane, (d): Radial cross section of panel (b). (e): Radial cross section of panel (c). 



2. Models. In this section, we introduce several neural field models used in the 
paper: our starting points are the models introduced by Laing et al. [42,43], in which 
the cortical space Vt is assumed to be the Euclidean plane M? . 

2.1. Integral model. The first model that we will consider is obtained from 
Eq. (1.1) assuming a translation- invariant, radially-symmetric kernel 



dt 



u{r,t) = —u{r 



.t)^ j w(\\v-v'\\)S{u{v',t))dv'^g{v) 



with sigmoidal firing rate (shown in Fig. 2.1(a)) 



radial connectivity function 

w ir) 



-br 



(6sinr + cosr). 



and external inhomogeneous input 

f ax^ + f3y^ 
g{r) = Go exp 



/i,6>> 0, 



|r||, b>0 



Go.aeR, a,/3>0. 



(2.1) 



(2.2) 



(2.3) 



(2.4) 



In the visual cortex regions of activation have been shown to have a Gaussian spread 
for radially-symmetric visual inputs [17], hence our choice for the function ^(r). The 
connectivity function is plotted in Fig. 2.1(b) on the Euclidean plane and as a ra- 
dial cross section in Fig. 2.1(d). An analytical expression for the Fourier transform 
of w cannot be obtained; we show the Fourier transform of Eq. (2.3) as computed 
numerically in the k = (/c^;, /c^y) -plane in Fig. 2.1(c) and as a radial cross section 
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in Fig. 2.1(e), where k = ||k||. The connectivity function (2.3) was proposed in 
models of working memory as a description of synaptic connections in the prefrontal 
cortex [34,44]. The connectivity describes local excitation and longer-range connec- 
tions that alternate between inhibition {w{r) < 0) and excitation {w{r) > 0). We 
argue that this type of connectivity pattern is also relevant to the study of patterns 
of activity in early visual areas like VI where there is a characteristic length scale 
associated with the average orientation hypercolumn width. It has been shown in 
anatomical studies that the number of lateral connections decay with distance, that 
the number of excitatory connections peak each hypercolumn width and the number 
of inhibitory connections peak each half- hypercolumn width [13]. The net effect is 
alternating bands of inhibition and excitation that decay with distance. This is also 
consistent with auto-correlations computed for the orientation selectivity map [49] 
given that connections tend to be reinforced between neurons with similar orientation 
preference. Henceforth the model (2.1)-(2.3) will be referred to as the integral model 
(IM). 

2.2. Fourth-order PDE approximation. In the cases when the Fourier trans- 
form of the synaptic kernel is a rational function, it is possible to derive an equivalent 
PDE formulation of Eq. (2.1) [42]. For simplicity, we will consider models without an 
external input ^'(r). If w{k) = P{k)/Q{k) with P and Q even functions in k where 
/c = ||k|| for k G M^, then the Fourier transform of Eq. (2.1) gives 



Q(fc) 



d 

—u{k,t) + u{k,t) 



= P{k){Sou){k,t). 



An inverse Fourier transform of the equation above leads to the desired PDE 

d 



^Q 



g^u{r,t) + u{r,t) 



= CpS[u{r,t)], 



where Cp and Cq are linear operators containing spatial derivatives of even order. 

Since the Fourier transform of the connectivity function (2.3) does not have an 
analytic expression, the integral model does not admit an equivalent PDE. However, 
we can approximate w with a function whose Fourier transform is rational and then 
derive an approximate PDE for the integral model. We specify the approximate 
connectivity function W4{r) through its Hankel Transform 

1 /^ ^ 

u)4{r) = t;- / swl{s)Jo{rs) ds, 
^TT Jo 

where 2^(/c) ^ w{k) is given by 

^(^) = B + ik^-Mr ^'-^^ 

and the coefficients A, 5, M are determined using a least-squares best fit algorithm 
(see Sec. 3.4.3 for further details). 

We compare the approximation with the original connectivity function in the 
physical and Fourier domains in Fig. 2.2(a) and (b). In physical space the two func- 
tions appear to be similar. The key qualitative difference is that in Fourier space 
2^(0) > 0, which is not consistent with the IM connectivity function, for which 
w{0) < 0. Biologically this means that w^ represents a globally excitatory connec- 
tivity function, whereas w represents a globally inhibitory connectivity function. We 
will see in Sec. 4.2.2 that it is necessary to increase the order of the polynomials in 
the numerator and denominator of Eq. (2.5) in order to accurately capture the sign 
at A: = 0. 

Starting from the expression for 2^, we derive the corresponding PDE, containing 
spatial derivatives up to the fourth order 



[5 + (M + A)2] 



— 2i(r,t) +2x(r,t) 



AS{u{r,t)), (2.6) 
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Fig. 2.2. Comparison of the PDE connectivity functions with w{r) given by Eq. (2.3) and its 
Fourier transform w(r) as computed numerically. Panels (a) and (h) show the Ath-order approxima- 
tion defined by Eq. (2.5) in the real and Fourier domains, respectively. Similarly, for the 8th-order 
approximation given by Eq. (2.7) in panels (c) and (d). For reference, we indicate the L2-norm of 
w, w — W4, w — wg and of their Fourier transforms. Parameters given at the beginning of Sec. 4-^- 



where the sigmoidal firing rate function is identical to the integral model case (2.2). 
In Eq. (2.6) we have denoted by A the standard Euclidean Laplacian. Henceforth, 
this model will be referred to as PDE4. 

2.3. Eight-order PDE approximation. In order to obtain a more accurate 
representation of the connectivity function (2.3), we repeat the steps outlined in 
Sec. 2.2 with the following approximation 



wi{k) 



-A{e-c){e-D) 

B^{k^- MY ' 



(2.7) 



where the values of A, 5, C, D and M are determined using a least-squares best 
fit algorithm; see Sec. 3.4.3 for further details. We compare this approximation with 
the original connectivity function in the physical and Fourier domains in Fig. 2.2(c) 
and (d). With the higher-order polynomials used here the approximation is more 
accurate and C^(0) > 0, consistent with IM (compare the 8th-order approximation as 
shown in Fig. 2.2(c) and (d) with the 4th-order approximation as shown in Fig. 2.2(a) 
and (b)). The synaptic function w^ leads to the following PDE, containing spatial 
derivatives up to the eight order 



[5 + (M + A)^] 



u{Y,t) -^u{r,t) 



= -[{A^C){A^D)]S{u{r,t)), (2.; 



where the sigmoidal firing rate function is identical to the integral model case. Hence- 
forth this model will be denoted as PDE8. 

3. Numerical methods. In this section, we review the numerical methods em- 
ployed for the computation of localised states in IM, PDE4 and PDE8. 

3.1. Integral model. Numerical computations of the IM (2.1)-(2.3) are per- 
formed discretizing a large but finite domain Q = [— L,L]^ with N evenly-distributed 
grid points in each spatial direction and imposing periodic boundary conditions. We 
approximate li on a grid Qn = { {^i^Vj) }^ -^i and collect the corresponding approx- 
imate values of li in a vector u 



Uij ^u{xi,yj), 



{ Uij }i 



N 



1>N^ 
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Similarly, we form vectors w, S(u), g G M^ for the approximations to w, S{u) and 
g^ respectively. Further, we introduce the discrete convolution, 

{u * v),j ^ T-^ [T(u) T(y)) (x„ %•), u * v = { (^/ * v),, }l^^^ G R^', (3.1) 

where we have denoted by T and T~^ the 2D Fourier Transform and its inverse, 
respectively. In summary, the discrete version of the evolution equation (2.1) is given 
by 

u = -u + w * S(u) + g. (3.2) 

This type of discretization has been applied before in direct numerical simulations 
of neural models (see, for instance [23,32]) even though it has not been used for 
numerical continuation. For smooth firing rate functions, the right-hand side can 
be evaluated accurately and efficiently using a Fast Fourier Transform (FFT) and 
its inverse (IFFT). In passing, we note that since the FFT of w can be performed 
and stored at the beginning of the computation, one function evaluation of the right- 
hand side requires just one FFT and one IFFT. Furthermore, standard de-aliasing 
techniques can be applied to the convolution operator if required [14]. 

Once a stable steady-state of Eq. (3.2) is found via direct numerical simulation, 
it is possible to continue it in one of the control parameters using standard numeri- 
cal continuation techniques. In previous studies of neural field equations, numerical 
continuation was performed on an equivalent PDF formulation of the integral system. 
A key observation is that path following can be applied directly to IM (or to similar 
models), employing FFTs and Newton-Krylov solvers [38,39]. Such methods do not 
require the formation of a Jacobian matrix, but rely only on Jacobian- vector multipli- 
cations: for IM, this is conveniently done using just a single application of FFT and 
IFFT. We remark that Newton-Krylov methods are often used in conjunction with 
sparse systems, but the performance of FFTs and IFFTs makes them a favourable 
choice for IM, even though the system is full. 

For numerical continuation of steady states of IM, we solve the system of algebraic 
equations 

F(u) = -u + w * S(u) + g = 0, (3.3) 

whose associated Jacobian- vector product is given by 

J(u)v = -v + w* (S'(u)v), u,vgM^', (3.4) 

where S'(u) = ^\^g{S' {u^x). • • • .S' {u^^)) G R^'><^'. We solve the system (3.3) 
using a Newton-GMRES solver implemented in Matlab and continue the solution 
with a secant method. Eigenvalue computations can also be performed using the 
Jacobian- vector products (3.4). Details of the numerical implementation and of the 
numerical parameters can be found in Sec. 3.4.1. 

Remark 3.1. The external input g guarantees that the system of algebraic equa- 
tions (3.3) is not translation invariant, even when the problem is complemented with 
periodic boundary conditions. Unless otherwise stated, we will use a negligible external 
input for the IM, so that Newton iterations can be applied directly to Eq. (3.3). We 
point out that a similar result can be obtained without imposing any external input, 
by perturbing the right-hand side with a term containing one extra unknown and then 
closing the system with a suitable phase condition [3, 15,46]. 

3.2. Fourth-order PDE approximation. In order to continue stationary lo- 
calised solutions to PDE4 with D6 symmetry, we use polar coordinates and pose a 
boundary- value problem on the sector ^i/q = { (r, ^) G M^ | < r < i^, < ^ < 7r/3 } 
with Neumann boundary conditions 

0=[B + {M + Af]u-AS{u), {r,e)enye 

0=Vu-n {r,0) e 3^1/6 ' 
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where the Laplacian operator A is expressed in polar coordinates. We discretise 
the system above using finite differences in r and a Fourier cohocation method in 0^ 
with Nr and Nq evenly spaced points, respectively, leading to a system of nonlinear 
algebraic equations of the form 

0= [m + (MI + L)2]u-AS(u) uGR^^^^ I,LgM^^^^><^^^^ (3.6) 

where I is the identity matrix. The Laplacian matrix L is formed explicitly, starting 
from differentiation matrices D^, D^^, T>qq for spatial derivatives with Neumann 
boundary conditions and then combining them with Kronecker products [3,46,56] 

L = D^^ (g) I^ + (R"^D^) (g) I^ + R"^ (g) Dee^ 

where R = diag(l,r2, . . . ,^Ar^) and Iq is the No-hy-No identity matrix. For purely 
radial patterns, we adapt the boundary- value problem so as to contain only the radial 
direction r and impose Neumann boundary conditions. Numerical continuation of 
the system (3.6) is performed with a secant method. Further details on the numerical 
implementation can be found in Sec. 3.4.2. 

3.3. Eight-order PDE approximation. For localised solutions of PDE8, we 
follow a similar approach to the on used for PDE4. In order to avoid the discretisation 
of 8th-order differential operators, we recast PDE8 as a system of two 4th-order PDEs 
and seek solutions to the following boundary- value problem 

0=(M + A)2ii-^, (r,^) G 1^1/6, 

=(M + A)'^ -^Bu^ A(A + C)(A + D)S{u), (r, 6) e Qi/q, 

0=Vii-n, {r,0) edQi/e, 

0=Vv-n, {r,0) edQi/Q. 



(3.7) 



Again, the discretisation of this system uses finite-differences in r and Fourier spectral 
collocation in 0. The example above shows also that, as we approximate w more 
accurately, the order of the underlying PDE increases and its numerical continuation 
becomes more demanding. A more convenient approach would be to use directly 
the integral form for the model, with connectivity function ws and proceed with a 
Newton-GMRES, solver. In this way, the computational cost would be the same as 
for IM. 

3.4. Implementation and numerical parameters. 

3.4.1. Integral model. Time simulations are carried out using a standard 4th 
order Runge-Kutta method with fixed step size. At each time step, the right-hand 
side of Eq. (3.2) is evaluated four times using an Nvidia Graphic Processing Unit 
(Tesla C2070). To compute the Discrete Fourier Transform we use the CUFFT library 
provided by NviDiA as part of its CUBA framework [48] . This software library allows 
us to easily exploit the parallelism of a GPU to obtain a fast implementation. The 
vector u is kept in the GPU global memory throughout a time step in order to avoid 
memory transfers and it is updated in parallel once the four stages of the Runge- 
Kutta scheme have been computed. Transfers between CPU and GPU memory only 
occur when the result of a time step needs to be saved to a file. Time simulations 
use a grid of 10^ x 10^ points and a stepsize of 0.5. Computation of each time step 
requires approximately 0.1 seconds. 

Numerical continuation for IM is done in Matlab using a in-house secant continu- 
ation code which employs a Newton-GMRES method for the nonlinear solves. Unless 
otherwise stated, we used a grid of 2^^ x 2^^ points and fixed an absolute tolerance 
of 10"^ for the nonlinear iterations. The Newton-GMRES solver uses the MATLAB 
in-built function gmres without preconditioners and with parameters restart = 20, 
tol = 10~^ and maxit = 10. Initial guesses are obtained directly from the Runge- 
Kutta time stepper, interpolating with MATLAB 's interp2 function where neces- 
sary. Stability computations are performed with Arnoldi iterations via MATLAB 's 



10 J. RANKIN, D. AVITABILE, J. BALADRON, G. FAYE AND D.J.B. LLOYD 

eigs function, passing the Jacobian- vector product (3.4) and computing (with the 
default tolerance) the first 20 eigenvalues with the largest real part. Computations 
are performed on a MacPro with a 3 GHz Quad-Core Intel Xeon processor employing 
exclusively the CPU. 

3.4.2. PDE4 and PDE8 models. Numerical continuation for the PDE models 
have been carried out with a secant code similar to the one used for IM, but using 
MATLAB's in-built function f solve for the nonlinear iterations. Unless otherwise 
stated, we used 300 grid points in the radial directions and 20 in the angular direction. 
We use the Levenberg-Marquardt algorithm implemented in f solve and set TolFun = 
10~^. The sparse Jacobians of these problems are formed and passed directly to the 
solver. Initial guesses for the continuation have been obtained using the expressions 
given in Eqs. (A.l) and (A. 2). Computations are performed on a standard laptop on 
a single core. 

3.4.3. Least-squares data fitting. Before comparing the connectivity func- 
tions W4{r) and ws{r) with w{r)^ it was necessary to tune the parameters in the 
definitions of wl and ws. In each case we performed a nonlinear least-squares opti- 
mization of the parameters using the Isqcurvef it function in Matlab. For W4^ the 
objective was to minimize the L2-norm of the difference between wl and iu whilst 
varying the parameters A, B and M in Eq. (2.5), where w is computed numerically 
using the Hankel transform at 300 points. Similarly, for ws^ the L2-norm of the dif- 
ference between ws and w was minimised whilst varying the paramaters A, 5, M, C, 
D in Eq. (2.7). For reference, the L2-norms of w and w are given above the panels in 
Fig. 2.2; the norm of the difference between the two functions plotted in each panel is 
also given. We note the largest Fourier mode of the connectivity dictates the location 
of bifurcations in terms of the sigmoid parameters and /i. By minimising the dif- 
ference between the connectivity functions in Fourier space, we expect to find similar 
behaviour for each connectivity over the same parameter ranges. On the other hand, 
if one minimises the difference between the connectivity functions in physical space 
and the amplitudes of the largest Fourier modes are not matched, bifurcations occur 
in different parameter ranges in each model and a direct comparison cannot be made. 

4. Numerical results. 

4.1. Convergence of the Nev^ton-GMRES solver. Since Newton-GMRES 
methods with pseudospectral evaluation of the right-hand side have not been used 
before for integral neural field models, we report briefiy on our solver. To test con- 
vergence, we perturbed a localised steady state of IM to obtain an initial guess (panel 
(b) of Fig. 4.1) and converge back to the original solution (panel (a) of Fig. 4.1) using 
our Newton-GMRES solver. In panel (c) we plot the relative residuals of each itera- 
tion, showing that we achieve convergence within a few nonlinear iterations. Similar 
convergence plots (not shown) are obtained for the numerical continuation, albeit so- 
lutions in that case are achieved with fewer iterations, owing to the more accurate 
initial guess provided by the secant predictor scheme. The experiment is repeated 
for various values of N: the convergence diagrams are indistinguishable from the one 
reported in panel (c). The wall time for the numerical experiment scales linearly with 
the number of unknowns A/"^, as reported in panel (d). We remark that, even with- 
out using any GPU acceleration and without enforcing explicit parallelisation in the 
CPU, the Newton-GMRES solver finds a solution to a fuh problem with 1,048,576 
unknowns in less than 40 seconds. 

4.2. Snaking behaviour of radial and D6 patterns. We now turn to the 
numerical continuation of localised states in IM, PDE4 and PDE8. The continuation 
parameter is the steepness /i of the sigmoidal firing rate, whereas the other parameters 
are fixed as follows: for IM, we choose = 5.6, b = 0.40, L = 60, Go = 10~^, a = 1.0, 
/3 = 1.0, a = VTO; for the PDE4 model, we use 6 = 5.6, R = 60, Go = with fitting 
parameters for wl in Eq. (2.5) given hy A = 1.225, B = 0.1398, M = 1.2183; for the 
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Fig. 4.1. Convergence of the Newton-GMRES solver. A stationary localised solution u^ of IM 
is shown in panel (a). The solution is perturbed to obtain an initial guess uq = u^ -\- 0.8 sin a; cos y^, 
shown in panel (b), and then the Newton-GMRES solver is used to converge back to w* . Convergence 
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where Wj is the solution at the jth itera- 



tion (panel (c)). The experiment is repeated for several values of N without any significant change 
to the convergence diagram, whereas the wall time for the numerical experiment scales linearly with 
the number of unknowns N'^ , as reported in panel (d). Parameters of IM: 6 = 5.6, fi = 2.5, b = 0.40; 
L = 60, Go = 4.0, a = 1.0, (3 = 4.0, a = 12.0. 



PDE8 model we use again ^ = 5.6, i? = 60, Go = with fitting parameters for wg 
in Eq. (2.7) given hy A = 0.8510, B = 0.6626, M = 0.6653, C = 0.3 and D = 10. 
We remark that translation invariance is removed in the IM by the negligible external 
input Go while in PDE4 and PDE8 this is achieved by the boundary conditions of 
the problems (3.5) and (3.7), so we choose Go = 0. In the present section we focus on 
the no- (or equivalent negligible-) input case which should be well understood before 
the addition of an input. In Sec. 5.2 we will provide examples of the model behaviour 
with input and discuss the implications. 

The bifurcation points in the diagrams presented in this section will be labelled as 
follows. F represents a fold bifurcation and P a spatial-symmetry-breaking bifurcation 
from a radial state. Superscripts indicate the symmetry properties of the bifurcation, 
where R represents a bifurcation on a branch with radially symmetric solutions and 
D6 represents a bifurcation on a branch of D6-symmetric solutions The labels / and 
r in the subscripts for fold bifurcations indicate whether the fold occurs on the left or 
right of the snaking structure. The indices n in the subscripts indicate the ordering 
moving up the snaking structure. On radial branches n corresponds to the number 
of rings around a central spot solution, for example, on the branch between Fj^ and 
F^ there is one ring around a central spot. On D6-branches there are n(n + l)/2 x 6 
additional spots glued around a central spot, for example, on the branch between F^^ 
and F^^ the solution has a total of 7 spots. 

4.2.1. PDE4 results and comparison with IM. We discuss both radially- 
and D6-symmetric localised solutions of PDE4 with connectivity function defined via 
Eq. (2.5). Other solution branches with different symmetry properties do exist but 
these are only discussed for IM in Sec. 4.3. An unstable radial spot solution bifurcates 
from the trivial state ix = at a Turing instability with /i- value to the right of the 
range shown in the subsequent diagrams. It is this unstable radial spot branch that 
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Fig. 4.2. Snaking of radially- and 1)6 -symmetric localised solutions in terms of fi for PDE4- 
Radial branches are grey curves and D6 branches are black; stable segments are solid and unstable 
segments are dashed. Bifurcations labelled F and P are discussed in the text, (a): Detail of radial 
branch from (b); points labelled al and a2 correspond to planar plots showing a stable spot solution 
and an unstable spot-with-ring solution^ respectively, (b): Global structure showing radial and D6 
branches; points labelled bl, b2 and b3 correspond to planar plots showing 7 -spot, IQ-spot and 37-spot 
solutions, respectively. Thin vertical lines discussed in text. Parameters given at the beginning of 
Sec. 4.2. 



appears in the bottom right of Fig. 4.2(a) and (b). 

Figure 4.2(b) shows the snaking structure for radial and D6 branches. We first 
focus on the radial branch, detail of which is shown in panel (a). A radial spot 
branch enters the diagram in the bottom-right-hand corner and undergoes a fold at 
FiQ. The radial spot solution existing on the branch segment between F^q and F^ 
is plotted in panel (al). This solution is stable between Fj^ and P^^. At P^^ on 
the radial branch a D6 instability produces the bifurcating branch of D6-symmetric 
solutions that leaves panel (a) in the bottom-left-hand corner. Beyond P^^ the radial 
branch is unstable and undergoes a further fold F^. After another fold Fj^ a ring has 
formed around the radial spot. The spot with ring solution existing on the branch 
segment between P^^ and P/^ is shown in panel (a2). The branch remains unstable 
and undergoes a series of further folds (P^f , P/27 ^^f ^ etc) adding additional rings as 
is shown in panel (b). 

The D6-symmetric branch that bifurcates from the radial branch at P^^ also 
undergoes a series of fold bifurcations as shown in Fig. 4.2(b). In this case a series 
of additional spots are added to the central spot in a configuration that preserves 
the D6-symmetry. Planar plots in panels (bl), (b2) and (b3) show the stable 7- 
spot, 19-spot and 37-spot solutions that exist on the branch segments between the 
fold pairs (P^?^P;?^), (P^?^PJ§^) and (P^?^P^^^), respectively. There are further 
intermediate stable branch segments between pairs of fold bifurcations that have not 
been labelled. On the stable segment that can be found between F^^ and P^^^, there 
exists a 13-spot solution for which one spot has been glued on the long edge of each 
of the six sides of the solution shown in panel (bl). Similarly, on the stable segment 
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Fig. 4.3. Comparison between the solution branches of PDE4 (greyscale) and IM (colour). 
Black curves and grey correspond to the PDE approximation and have the line style and labelling 
conventions as Fig. 4. 2. Corresponding radial and D6 solution branches from IM are plotted as 
dash-dot curves in red and blue, respectively; stability is not indicated and folds not labelled on these 
branches, (a): Detail close to the spatial- symmetry-breaking bifurcation P^^ ; the corresponding 
bifurcation for IM is labelled Pj^ • (b): Detail of the D6 snaking structure. Parameters given at 
the beginning of Sec. 4--2- 



that can be found between F^^ and Fj^^ ^ there exists a 25-spot solution for which 
two spots have been glued on the long edge of each of the six sides of the solution 
shown in panel (b2). 

The same radially- and D6-symmetric branches shown in the previous section 
have been computed for IM. Here we test the accuracy of PDE4 both qualitatively 
in terms of the types of solutions produced and their bifurcations, and quantitatively 
in terms of the parameter ranges for which the different solution types persist. We 
are also interested to see whether the relative ranges of existence for different types 
of solution is consistent between PDE4 and IM. 

Figure 4.3(a) and (b) both show detail from Fig. 4.2(b) with the same curves 
reproduced for PDE4 with the same line style and labelling conventions. Also plotted 
(in colour) are the equivalent curves computed for IM, where the equivalent of P^^ 
in IM is P^ . The first major point to make is that in terms of the types of solution 
encountered, the series of bifurcations encountered and the stability of each branch 
segment, there is an exact agreement between PDE4 and IM. Furthermore, the quan- 
titative agreement on the radial branch is good up until the fold point F^ . Above F^ , 
the radial branch for PDE4 makes a large excursion away from the IM branch and the 
branches remain well separated as the snaking continues; see panel (a). Similarly, for 
the D6 branch the level of agreement is good up to F^^ above which PDE4 branch 
deviates and remains well separated from the IM branch; see panel (b). The range 
of existence in ji for each stable branch segment on the D6 branch is greatly under 
estimated by PDE4. 

We now highlight a key qualitative difference between the bifurcation diagrams 
for PDE4 and IM. In IM, there is a range of /i G [2.5, 3.0] for which the stable branch 
segments corresponding to 7-spot, 19-spot and 37-spot all overlap. This is not the 
case for PDE4, in particular, the branch segments corresponding to stable 7-spot and 
37-spot solutions between the fold-pairs (F^^.F^^) and (F^^.F^^) do not overlap. 
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Fig. 4.4. Snaking of radially- and 'D6-synimetric localised solutions in terms of fi for PDE8 
(grey scale) and comparison with IM (colour). Black and grey curves correspond to the PDE approx- 
imation and have the line style and labelling conventions as Fig. Jf..2. Corresponding radial and D6 
solution branches from the full model are plotted as dash-dot curves in red and blue, respectively; 
stability is not indicated and folds not labelled on these branches, (a): Global snaking diagram for 
PDE8. (b): Detail close to the spatial- symmetry-breaking bifurcation P^^ ; the corresponding bifur- 
cation in the full model is labelled Pj^ ■ (c): Detail of the D6 snaking structure. Parameters given 
at the beginning of Sec. 4-2- 



This can be seen by the fact that F^^ occurs at a smaher /i-value (indicated by 
the first thin vertical hne in Fig. 4.2(b)) than F^^ (indicated by the second thin 
vertical line). This organisation of the solutions in parameter space is qualitatively 
inconsistent with IM. 

4.2.2. PDE8 results and comparison with IM. Figure 4.4(a) shows branches 
of both radially- and D 6- symmetric solutions of PDE8. Globally the bifurcation di- 
agram is the same as that of PDE4 in terms of the types of solution observed and 
the sequence of bifurcation encountered. There is an important difference between 
PDE4 and PDE8 in terms of the organisation in parameter space of the solution 
branches. For PDE4, there is no overlap in parameter ranges for which the 7- and 
37-spot branches are stable; see the two vertical lines in Fig. 4.2 and note that F^^ 
occurs before Ff^^ in this case. In the PDE8 case, as shown in Fig. 4.4(a), there 
is an overlap in the parameter ranges as indicated by the grey shaded region. This 
organisation of the solution branches in parameter space is now consistent with the 
full model as shown in panel (c). Indeed PDE8 provides better agreement with IM; 
the branches remain close for both the radial and the D6 branches as we move up the 
snaking structure as can be seen in panels (b) and (c). We note that when compared 
with IM, the upper section of the radial branch occurs at smaller values of /i for PDE4 
and at larger values of /i for PDE8; compare Fig. 4.2(a) with Fig. 4.4(b). 

4.3. Snaking of D2, D3 and D4 patterns in IM. In this section we discuss 
several patterns that possess neither radial nor D6 symmetry. For non-radial pat- 
terns discussed so far in this article, see Fig. 4.2(bl)-(b3), the individual spots lie on 
a regular hexagonal lattice. However, there exist an infinite number of stable config- 
urations that conform to the same lattice spacing but without the full D6 symmetry; 
here we present two such examples. We also present one further example of a stable 
configuration that does not conform to the regular hexagonal lattice. 

4.3.1. Three-spot pattern (D3). Figure 4.5 shows the snaking of D3-symmetric 
patterns about a three-spot solution. The unstable branch that enters panel (a) in 
the bottom-right-hand corner reconnects to the trivial state ix = at the Turing in- 
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Fig. 4.5. Snaking of patterns with D3 symmetry built around a triangular- organised, three- 
spot solution. Line style conventions for stability as in Fig. 4- 2- Parameters for IM as given at the 
beginning of Sec. 4. 2. 
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Fig. 4.6. Snaking of patterns with D2 symmetry built around a two-spot solution. Line style 
conventions for stability as in Fig. 4. 2. Parameters for IM given at the beginning of Sec. 4. 2. 



stability (not shown). The panels (al)-(a5) show stable solutions on the first five 
full excursions in \i of the snaking structure. The existence of three-spot (see panel 
(al)) and twelve-spot (see panel (a3)) solutions for PDE4 with connectivity given by 
Eq. (2.5) was shown in [43]. Here we have shown that these solutions exist in IM and 
that they form part of a larger snaking structure. 

4.3.2. Two-spot pattern (D2). Similarly, there is an unstable two-spot solu- 
tion that connects to the trivial state i^ = at the Turing instability (not shown). 
Figure 4.6(a) shows that this solution also undergoes a sequence of fold bifurcations 
giving rise to larger D2-symmetric patterns. We note that the spacing between the 
spots in these patterns still conforms to the regular hexagonal lattice. The pan- 
els (al)-(a5) show solutions on the first five stable branch segments moving up the 
snaking structure; we note that the pattern (a3) is on an intermediate branch that 
does not make a full excursion in /i. 

4.3.3. Quincunx pattern (D4). We show in Fig. 4.7(a2) a stable configuration 
that lies on a square lattice but with a spacing between the spot peaks that is double 
that of the hexagonal-lattice solutions encountered thus far. The solution is formed 
by five spots that interact at the first excitatory peak away from in the connectivity 
function (2.3) as shown in Fig. 2.1(d). The configuration of four spots forming a square 
with an additional spot in the center is typically referred to as a quincunx pattern that 
is found, for example, on dice and dominoes. As shown in Fig. 4.6(a) these solutions 
exist on an isola in parameter space where other branches, see panels (al) and (a3), 
are unstable. Although the pattern does not undergo snaking-type behaviour it may 
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Fig. 4.7. Isola of five- spot solution organised in a quincunx configuration. Line style conven- 
tions for stability as in Fig. 4-^- Parameters for IM given at the beginning of Sec. Jf.-^- 



be possible to construct larger patterns on the double-spaced square lattice. 
5. Discussion. 

5.1. Summary. This paper explores patterns of localised activity in the neural 
field equation posed on the Euclidean plane with a smooth firing rate function. The 
choice of connectivity function is an important factor in determining whether, the 
localised behaviour found is restricted to individual spots, or whether multiple inter- 
acting spots can form coherent localised patterns. In [30] localised states were studied 
in a model with a radially-symmetric wizard-hat connectivity function describing local 
excitation and lateral inhibitions in the Euclidean plane. When spot solutions were 
tracked using numerical continuation no snaking behaviour was observed, i.e. the 
only steady states found consisted of a single spot. The radially symmetric connec- 
tivity function studied here and shown in Fig. 2.1(d)) features local excitation, lateral 
inhibition and long-range bands of excitation that decay with distance. Indeed, the 
distance between excitation peaks fixes a spatial scale that allows for regular spatial 
interactions and the formation of larger patterns of activity. In [43] it was shown that 
multiple-spot patterns could be obtained all with the common property of the peaks 
lying on a regular hexagonal lattice with spacing determined by the connectivity func- 
tion. One of the main aims in the present article was to show how these solutions are 
connected in parameter space and how patterns with varying spatial extent grow via 
the mechanism of homoclinic snaking. 

The results presented in [43] relied on working with an approximated connectivity 
function (see Fig. 2.2(a) and (b)) that allowed for solutions the full integral neural 
field equation to be studied in an equivalent fourth-order PDF . The initial parts of 
the results section are concerned with the relative agreement between solutions to 
the integral model and equivalent PDF formulations with approximated connectivity 
functions. We pursue the problem numerically by investigating the level of agreement 
in terms of an entire bifurcation diagram rather than comparing individual solutions 
at fixed parameter values. We compared both radially symmetric and D6-symmetric 
solution branches and found that the qualitative difference of the zero-mode in the 
Fourier domain for the approximated connectivity used in the fourth-order model 
(see Fig. 2.2(b)) led to significant discrepancies in the existence ranges of solution 
branches, in particular for solutions with a larger spatial extent. We demonstrated 
that the improved approximation of the connectivity function shown in Fig. 2.2(c) 
and (d), leading to an eighth-order PDF, provides a better agreement across the full 
bifurcation diagram. In particular, the eigth-order model captures the key feature of 
there being a specific parameter range in which multiple solutions coexist, each solu- 
tion with a different spatial extent. It was not possible to capture this feature with 
the fourth-order approximation. We conclude that, although equivalent PDF formu- 
lations have proved to be a useful tool for the study of neural fields, it is important 
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Fig. 5.1. Model simulations with input given by Eq. (2.4) with a = (3 = 1. The simulations are 
performed with /j, = 2.4. (a) The Gaussian-bump input plotted with cr = 10 and Gq = 6 (Gq = 1.5 
used in the simulations) . (b): Case with a = 9.0, the model converges to a 7s-spot solution that is 
a modified version of the state shown in Fig. 4-2(bl) (c): Case with a = 9.5, the model converges 
to a 12s-spot solutions that is a modified version of the state shown in Fig. 4-5(a3) (d): Case with 
(7 = 10.0, the model converges to a lAs-spot solutions that is a modified version of the state shown 
in Fig. 4-6(a5). Other model parameters (for IM) given at the beginning of Sec. 4-2. 



to ensure close agreement between the connectivity functions in the Fourier domain. 
Increasing the order of the PDEs used ahows for improvements in this agreement. We 
beheve that, while converting the integral formulation to higher-order PDEs could 
be useful in analytical studies, numerical calculations of these systems should be ap- 
proached without resorting to PDE formulation where possible. In passing, we point 
out that the methodology proposed here for the integral model is applicable to inho- 
mogeneous synaptic kernels, provided that the convolution structure of the integral is 
preserved. Furthermore, we point out that here we have used the standard Newton- 
GMRES method mainly for its simplicity, but more sophisticated choices are also 
possible [39]. 

Having investigated radial and D6 solutions with PDE formulations and com- 
pared the results with the full integral model, we proceeded to give an account of 
other types of solutions that, when path- followed with numerical continuation lead 
to patterns with different underlying symmetry properties. We worked with the full 
integral model and showed that patterns with D2 and D3 symmetry also give rise to 
snaking behaviour, generating spatial patterns with variable spatial extent. All the 
solutions of this type, including those shown earlier with D6-symmetry, have the com- 
mon feature of the individual spots lying on a regular hexagonal lattice. Furthermore, 
these solutions of variable spatial extent exist within roughly the same ranges of the 
parameter ji. As the snaking diagram is ascended, new spots are glued to long edges 
of the pattern in a regular fashion. We note the existence of an arbitrary number 
of intermediate branches not shown in our bifurcation diagrams where symmetries 
can be broken via the (simultaneous) addition or subtraction of one (or more spots). 
We expect these solutions to lie on intermediate branches that are stable over ranges 
smaller than the full excursions in the main snaking structure. 

5.2. Localised patterns with input. The numerical bifurcation analysis pre- 
sented in this paper opens up the possibility to investigate the spread of cortical 
activity with inputs. An important principle in studying the neural field equation is 
that weak inputs to the equations should drive the system to states that are already 
solutions to the underlying equations without input. The bifurcation study presented 
here allows us to identify the types of solution that we may expect to encounter for 
localised inputs and the relevant parameter regime in which they occur. Two criteria 
need to be satisfied when identifying a suitable operating regime, 1) the model should 
only produce the trivial homogeneous state ii = before an input is introduced and 
2) when an input is introduced, the model should be driven to one of the underly- 
ing non-trivial solutions. We have shown that, for the full integral model, there is 
an accumulation of fold bifurcations at around ji = 2A representing the first point 
for which localised patterns can be observed. Both criteria are satisfied when the 
model is operated just before these fold points. Introducing a weak input the sys- 
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Fig. 5.2. Symmetry breaking of weakly unstable solution on a square lattice at jjl = 3.2. (a): 
Evolution of L2 norm, (b): Initial condition given by Eq. (A. 3). (c): Weakly unstable solution, 
(d): Stable solution. Parameters (for IM) given at the beginning of Sec. J^..2. 



tern can be driven to states that have a spatial extent corresponding to that of the 
input. Figure 5.1(a) shows the profile of the input used in a series of simulations 
that were initiated with the ix = state plus a small random perturbation across 
the entire spatial domain. Figures 5.1(b), (c) and (d) show the result of three 150 
time unit simulations, each with different spatial extent for the input. In each case, 
the trivial state ix = is no longer stable and the system naturally selects one of 
the solutions described in the early bifurcation analysis. As the spatial extent of the 
input increases, the size of the pattern selected by the model increases and this is an 
important consequence of the corresponding solutions existing as part of a snaking 
structure. The computational framework presented in this article will allow for the 
relationship between model inputs and spatially localised patterns to be investigated 
in future work. 



5.3. Patterns on other lattices and parametric forcing. The multi-spot 
solutions described in this article all have the common feature of the activated peaks 
falling onto a regular hexagonal lattice. Indeed, this has been found to be the default 
way for the radial symmetry to be broken in pattern forming systems, notably the 
archetypal Swift-Hohenberg equation [46]. We also investigated whether it is pos- 
sible to find other states that do not conform to the regular hexagonal lattice. In 
Fig. 4.7 an example of a solution with D4 symmetry was shown that consists of five 
spots interacting at double the standard separation between excitation peaks. We 
found this solution to exist on an isola and not undergo snaking so as to find larger 
patterns tiling the plane with the same spacing. We also attempted to converge so- 
lutions with D4 symmetry that have the regular spacing between peaks. A suitable 
initial condition to find such solutions is given by Eq. (A. 3) chosen such that there 
is a depression at (x^y) = (0, 0) and the surrounding square-lattice pattern decays 
away from the origin, see Fig. 5.2(b). In an appropriate parameter range these states 
appear to converge to stable patterns. However, we found the patterns to be weakly 
unstable, finally converging to a pattern on a hexagonal lattice after a long transient. 
Figure 5.2(a) shows a time-course of the L2 norm from a simulation with the initial 
condition shown in panel (b). The model reaches the apparently stable D4 config- 
uration after approximately 30 time units (panel (c)) before finally converging after 
230 time units to the D2-symmetric pattern (panel (d)) that was previously identi- 
fied in Fig. 4.6(a2). It would be possible to stabilise such weakly unstable solutions 
with the use of parametric forcing, by introducing small modulations that encourage 
interactions on a fixed lattice in cortical space. In the neural field equations this 
is typically referred to as inhomogeneous neural media; travelling waves, travelling 
fronts, periodic patterns and pulsating fronts have been studied with such modula- 
tions [7,20,24]. The study of localised states in this context would be an interesting 
future direction. Furthermore, stable quasi-periodic patterns have been obtained in 
the Swift-Hohenberg equation through parametric forcing [37]. The question of intro- 
ducing orientation-preference tilings on square and hexagonal lattices was addressed 
in [4]. However, one is restricted in the number of orientations that can be equally 
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represented with such tihngs, (two and three, respectively). Parametric forcing on a 
quasi-periodic lattice is of particular interest in the neural field equations as this could 
allow for near-continuous representations of features in a model without an abstracted 
feature space. 

6. Conclusions. The organisation in parameter space of localised structures 
consisting of multiple spots has been revealed for the first time in planar neural field 
equations. In order to find such behaviour, one must choose a connectivity function 
with an excitatory peak away from the origin that fixes a regular spatial scale of inter- 
actions between spots. As localised solutions are path- followed using numerical con- 
tinuation we find that these structures grow in a series of fold bifurcations through the 
mechanism of homoclinic snaking that has been well-studied in the Swift-Hohenberg 
equation. A numerical strategy has been proposed to perform a numerical bifurcation 
analysis without resorting to a PDE formulation, but taking advantage of matrix-free 
Newton-Krylov nonlinear solvers combined with a pseudospectral evaluation of the 
right-hand side. The novel application of these methods to the neural field equations 
allowed for numerical continuation to be applied to the full integral form of the model. 
Previous studies in 2D have relied exclusively on PDE approximations of the connec- 
tivity functions; here we demonstrated that these approximations can give a very 
close agreement with the full integral model if a sufficiently high-order approximation 
is taken. The numerical schemes presented here will allow for future studies of the 
neural field equations to use connectivity functions defined either directly in the real 
domain or the Fourier domain without recourse to PDE methods, provided that the 
sigmoidal firing rate be smooth and that the integral formulation can be expressed as 
a convolution (this extends also to inhomogeneous firing rates). 

The neural field studied in the present paper can be considered as a model of 
the visual cortex and the localised patterns studied without inputs can be related 
to visual hallucinations that can be localised in the visual field [54]. Furthermore, 
we have shown that the localised states computed in our bifurcation analysis are 
exactly the types of solutions selected by the model in the presence of weak inputs. 
The persistence of these localised structures in the presence of a model input is new. 
This future direction will be of particular interest for the study of localised patterns 
of activity that have been observed in the primary visual cortex [17] with localised 
visual input. 
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Appendix A. Analytic expressions for initial conditions. Figure 1.2(a) 
shows the initial condition given by 

u{x,y) = Aexp f — j, (A.l) 

with A = 6 and L = 5.77; subsequent panels (b) and (c) show a transient state after 
1 time unit and the stable steady state after 15 time units. Figure 1.2(d) shows the 
initial condition given by 



u{x, y) = Aexp( — 



I V3 \ 1 V3 

COS(x) + COS -X + -T^y + COS — 7^^ + 



(A.2) 

with A = 2 and L = 100; subsequent panels (e) and (f) show a transient state after 
1 time unit and the stable steady state after 15 time units. Figure 5.2(b) shows the 
initial condition given by 

u{x ^ y) = 2 exjp I I (— cosx — sin?/), (^-3) 
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with A = 2 and L = 65; subsequent panels (c) and (d) show transient states after 
200 and 300 time units, respectively. The L2-norm is plotted over this time course in 
panel (a). 
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